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Abstract 

Several approaches are presented to identify an experimental system model 
directly from frequency response data. The formulation uses a matrix-fraction 
description as the model structure. Frequency weighting such as exponential 
weighting is introduced to solve a weighted least-squares problem to obtain 
the coefficient matrices for the matrix-fraction description. A multi-variable 
state-space model can then be formed using the coefficient matrices of the 
matrix-fraction description. Three different approaches are introduced to fine- 
tune the model using nonlinear programming methods to mimmize the desired 
cost function. The first method uses an eigenvalue assignment technique to 
reassign a subset of system poles to improve the identified model. The second 
method deals with the model in the real Schur or modal form, reassigns a 
subset of system poles, and adjusts the columns (rows) of the input (output) 
influence matrix using a nonlinear optimizer. The third method also optimizes 
a subset of poles, but the input and output influence matrices are refined at 
every optimization step through least-squares procedures. 


1 Introduction 


One major objective of system identification is to provide mathematical models for 
dynamics and control analysis and designs. However, models of systems can have 
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various forms, such as transfer functions, differential or difference equations, and 
state-space equations. A frequency-domain state-space identification method [1-5] 
provides a state-space model of a linear system from frequency response data. 

The method called the State-Space Frequency Domain (SSFD) identification al- 
gorithm [2] can estimate Markov parameters (pulse respond) from the frequency 
response function (FRF) without window distortion when an arbitrary frequency 
weighting is used to shape the estimation error. The method uses a rational matrix 
fraction description (the ratio of a matrix polynomial and a monic scalar polynomial 
enominator) to curve-fit the frequency data and compute the Markov parameters 
rom FRF. The curvetting problem must be solved either by nonlinear optimization 
techniques or by linear approximate algorithms with several iterations. To obtain the 
state-space models from the Markov parameters, the Eigensystem Realization Algo- 
rithm (ERA) or its variant ERA/DC is used [5], 

Frequency domain methods presented in Refs. [3, 4, 5] start with identifying a 
eft matrix-fraction description (LMFD) of the transfer function matrix. The ad- 
vantage of using the LMFD, as an intermediate model between the data and the 
esire final state-space model, is that from frequency response data to the LMFD 
is a linear least-squares problem, which is easy to solve. This method works quite 
well when the frequency response data are fairly accurate; however, it might yield 
unstable, erroneous models if the data contains too much distortion and/or error 
ata distortion in the frequency domain is caused by a number of factors; limited 
sampling frequency, filters to remove noise, and lack of periodicity. This data dis- 
tortion often cause unstable modes to be present in the identified system model An 
improved method was introduced in Ref. [6] to deal with the problem when data 
istortion is present. The idea is to stabilize or remove the unstable modes before 
expanding the matrix-fraction description (MFD) into the Markov parameters (pulse 

responses). This approach avoids introducing unstable modes while still maintaining 
the frequency response close to the data. 

In this paper, exponential frequency weighting [2, 7] is used to solve a weighted 
least-squares problem for the LMFD coefficient matrices. A multi-variable state-space 
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model is then realized from the LMFD coefficient matrices. To improve the identified 
model, nonlinear programming methods [8] are used to fine-tune the model param- 
eters. There are three different formulations introduced in this paper for parameter 
optimization. In all three formulations, the objective function is defined as the error 
between the actual FRF and the synthesized FRF using the identified space-space 
model. The first formulation uses a general system realization, and utilizes nonlinear 
programming along with an eigenvalue assignment [9 — 11] technique to optimize a 
subset of system poles. The second formulation deals with system realizations in 
the real Schur or modal forms, and uses a subset of system poles, as well as some 
coefficients to adjust the columns (rows) of the input (output) influence matrix for 
parameter optimization. The third formulation is similar to the second, but the input 
and output influence matrices are computed at every optimization step through least- 
squares procedures. Experimental data from a NASA testbed with fifteen inputs and 
fourteen outputs arc used with a total of two hundreds and ten transfer functions to 
demonstrate the concepts proposed in this paper. 

2 Weighted Least- Squares Method 

Given the system frequency response function G(z k ) at the frequency point z k , con- 
sider the left matrix-fraction 

G(z k ) = a~ 1 (z k )f3(z k ) (1) 

where 

OciyZ k ) = I m + OL\Z k 1 H OL p Z k P (2) 

0{z k ) = Po + Pl z k 1 d +/^P' 2 fc P (3) 

are matrix polynomials with I m being an identity matrix of order m. Every a, is 

an m x m real square matrix and each $ is an m x r real rectangular matrix. The 

factorization in Eq. (1) is not unique. For convenience and simplicity, one can choose 
the orders of both polynomials to be equal to p. 
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Pre-multiplying Eq. (1) by a(z k ) produces 

a(z k )G(z k ) = f3(z k ) 

which can be rearranged into 


G {z k ) = -ct x G(z k )zZ x <* P G{z k )z k p 

+ Po + Piz k l 4 + P p z k p 

or 


( 4 ) 


(5) 


G(z k ) = ®g k (g) 

where the matrix 0, of dimension m x [p(m + r) + r], and the matrix Q k , of dimension 
\p{m 4- r) + r] x r, are defined as 

0 = ai • • • a p p 0 Pi fa . . . /J p j (7) 

G(z k )z k 1 ' 

G(z k )z k p 

I r ( 8 ) 

IrZ k l 
I r z k p 

Here, I r is an r x r identity matrix. With G(z k ) and z k l known, Eq. (5) or (6) is a 
linear equation. Because G(z k ) is known at z k = (k = 1 , . . . there are £ 

equations available. 

The parameter matrix 0 in Eq. (6) is a real matrix whereas G(z k ) and Q k are 
both complex matrices. Thus Eq. (6) is a complex matrix equation with a total of £ 
complex equations . Let us define 

G k = [Real(G(2*)) Imag(G( 2fc ))] and Q k = [Real(£ fc ) Imag(^)] (9) 
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Equation (6) may be rewritten as 


( 10 ) 


G k = eg k 

Equation (10) is a real matrix equation consisting of 2£ linear equations for computing 
the parameter matrix 0. The matrix G k at the frequency point k is an m x 2r matrix, 
whereas Gk is a [p{ m + r) + r] x 2r matrix. 


2.1 Recursive Formulation 

To solve Eq. (10), let us first define a weighted cost function to be minimized as 

J(0, k) = II ©(&){?<-» - Ge-i 111 • 

»=i 

where 0 < w < 1 is a forgetting factor weighting the frequency data. The data at 
the lowest frequency point is given unit weight, but data that is k frequency points 
higher is weighted by w k . The method is commonly called exponential forgetting. 
The cost function defined in Eq. (11) is motivated by the fact that accelerometers arc 
commonly used as the measurement device in structural testing. The corresponding 
frequency response functions have better response levels in the high frequency range. 
Identifying lower frequency information in the presence of measurement noise becomes 
a problem. One way to solve this problem is to weight more the lower frequency region. 
On the other hand, displacement sensors have better response capability for the low 
frequency region. For this case, the forgetting factor may be switched to weight the 
high frequency region more than the lower frequency region. The form of Eq. (11) is 
unchanged except for the index £ i is replaced by i. 

Using recursive least squares, the solution that minimizes Eq. (11) is 

G(k) = <j>(k)P(k) ( 12 ) 


where 


(p(k) 

P(k) 


1=1 

^w^Ge-iGh 

Lt=l 


i -1 


(13) 

(14) 
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The matrix P(k) is the inverse of the frequency data correlation matrix weighted by 
the forgetting factor w. Note that the matrix P(k) is positive definite. Prom the 
definition of P(k), application of the matrix inversion lemma yields 

m = 

A=1 


= ^ p (k - 1) - ~P(k - 1 )S t ^IC(k - 1) 

where 

K(k - 1 ) = [wi 2r + gj_ t p(k - 1 )9t- k ]- l gl k p(k - 1 ) 

Similarly, the quantity (p{k + 1) at frequency point k + 1 can be written as 

m = 

1=1 

= w<j>(k - 1) + G e _ k gJ_ k 

Substituting Eqs. (15) and (16) into Eq. (12) yield the parameter 0(/fc) 
O(k) = 4>(k)P(k) 


(15) 

(16) 

(17) 


M* - 1) + G e _ k gj_ k ] - 1) - ±P(k - 1 )Q t _ k K{k - l)" 

= Q(k - 1 ) + [G(_ k - e(k - 1 )g e -k] tC(k - 1 ) 

where the last equality results from the following 

6tk [ l P{k “ - i p ( k - 1 )S,- k IC(k - 1 )] 

= w + ^< -t ptk ~ !)&-*] - BL t P(k - lc(k - 1 ) 
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The parameter 0(fc) which minimizes the cost function is given recursively by 

G(- k = [Real(G(z*_ fc )) Imag(G(**_ fe ))] 

G(z e -k)z^ k 

G(z e - k )zif k 
I r 

Ir z l-k 
Ir z e-k 

Ge-k = [Real(&_ fc ) Imag(ft_ fe )] 

K{k - 1) = [whr + gJ_ k P(k - 1 )G t -kY'Gl- k P{k - 1) 

G e -k = 0(k-l)Ge-k 

pm = -P(k-l)--P(k-l)Gt-klC{k-l) 

U) w 

6(fc) = e (k - 1) + [Ge-k - Ge-k] IC(k - 1) 

At any specific frequency point k, the m x 2r FRF matrix Gf~k is given, where r is the 
number of inputs, m is the number of outputs, and i is the total number of frequency 
points to be used for the identification process. The [p(r + m) + r] x r complex matrix 
Gk is then computed where p is an integer large enough to satisfy the constraint 
pm > n (the system order). The 2 r x [p(r + m) + r] matrix lC(k - 1) is the update 
gain determined by the matrix P(k - 1) of dimension [p(r + m) + r] x \p(r + m) + r], 
the matrix Ge-k of dimension [p(r + m) + r] x 2 r, and the scalar w. The initial values 
of P(0) and 0(0) can be arbitrarily assigned, normally, P(0) and 0(0) are assigned 
as dl and 0, respectively, where d is a large positive number, I is an identity matrix 
dimensioned \p(r + m) + r] x [p(r + m) + r], and 0 is a zero matrix dimensioned 
m x [p(r + m) + r]. 

In the recursive process, special care must be taken to ensure that both matrices 
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[wl 2r + Gji k P(k - l)G?-kY l and P(k) must be symmetric and positive definite. In 
theory, inverting a positive-definite matrix results in another positive-definite matrix. 
In practice, any numerical error in the matrix inversion process at any step may 
accumulate large enough errors to destroy its symmetry and positive-definiteness. It 
eventually leads to an unstable solution, i.e., the parameter © will not converge to a 
constant value when k -» £. One simple way to eliminate the inversion problem is 
to take only the symmetric part of the inverted matrix at every recursion step, i.e., 

2 { [ U ’ I‘ 2 r Q P ( k — 1 )Gi~k\ + (\wI 2 r + Gj_ k P(k — l)cj^_fc] -1 ) T }. This will guarantee 
that the inverted matrix is symmetric. 

2.2 Batch Formulation 

Recursive approaches are better suited for computations in real time, i.e., parameters 
are computed as data becomes available. Often, experimental data from a completed 
test is available which allows all calculations to be performed at once. A batch version 
is presented in this section. Stacking up the 2£ equations in Eq. (10) yields 

6 = eg (20) 

where 

G = G 0 Gi * • ■ G e ] 

G = [ Go Gi • • • Q e ] (21) 

The cost function J shown in Eq. (11) is minimized by solving the least-squares 
solution for 0 according to Eqs. (12), (13), (14) with k = £, 

e = GGimr 1 ( 22 ) 

where 

Gw = Go wG\ • • • w e Ge ] (23) 

The subscript w associated with G w signifies the forgetting factor w inserted into Q 
with an appropriate power at each frequency point. 

The weighting w e for the highest frequency at the frequency point i can be quite 
small depending on the length £ of the data and the choice of the forgetting factor w. 
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Fox example, w e « 4.3 x 10' 5 with l = 1000 and w = 0.99. Unless the amplitudes of 
those frequencies near the highest frequency are in the order of 10~ 5 , their contribution 
to the identification process may become negligible. Using accelerometers, the ratio 
of the highest frequency to the lowest frequency can be as high as 10 3 to 10 5 . For 
this case, the forgetting factor used in Eq. (23) is indeed a good weighting technique 
to perform a better low-frequency identification. 

On the other hand, one may prefer to have freedom of choosing a weighting factor. 
A slight modification of Eq. (23) will provide such freedom, i.e., 


Go w\G\ 


wtGe 1 


(24) 


The quantities w u w 2 , ■ ■ ■ ,Wi, can be all independent. They may be randomly or 
specifically chosen. Some obvious choices include 

w k = e“ 10(1_fc) ^, w k = -, w k = £ 2 ? k = l,2,...,£ 

For the case where the low frequency resolution is better than the high frequency 

resolution, the weighting must be reversed. 

Substituting Eq. (24) in Eq. (22) and solving for the parameter 0 that minimizes 

the following cost function, 

J(0, £) = ^ w i II - Hi 


yields results similar to Eq. (11) except for the weighting factor. 

In the next section, optimization-based approaches to further improve the least- 

square solution are discussed. 


3 Nonlinear Optimization 

Another approach to enhance the identified model is to use nonlinear programming 
to tune the model parameters obtained from the solution to Eq. (10). Once the 
solution, represented by the parameter matrix 0, is computed using Eq. (18) or 
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Eq. (22), a state-space realization is determined. The state-space realization can be 
m any canonical form such as Schur form, modal form, Jordan form, observable form 
etc. There are three different formulations considered in this paper for parameter 
optimization. The first formulation deals with a general system realization, and 
utilizes nonlinear programming along with an eigenvalue assignment technique to 
optimize a subset of system poles to improve the agreement between the measured 
transfer function and the identified model. The second formulation deals with system 
realizations m the real Schur or modal forms, and uses a subset of system poles, as well 
as some newly defined parameters to adjust the columns (rows) of the input (output) 
influence matrix, as optimization parameters. The third formulation is similar to the 
second, but the input and output influence matrices are not directly adjusted by the 
optimizer, rather, they are computed at every optimization step through least-squares 
procedures. In the first formulation, all system poles are reassigned simultaneously 
e desired values given by the optimizer, via an eigenvalue assignment technique, 
n the second and third approaches , each pole is individually reassigned by the 


3.1 Parameter optimization: Eigenvalue Assignment 

In this formulation, a subset of system poles are used as optimization parameters 
to minimize a cost function, which measures the difference between the experimen- 
tal transfer function and the identified transfer function over frequency range of in- 
terest. Of course, the direct approach would be to use the elements of the state 
matrix directly, with equality constraints to reassign the poles. However, this ap- 
proach is computationally expensive since it requires too many design parameters. 
Let (A 0 ,B,C\D) represent an initial realization for the identified system. As men- 
tioned earlier, to determine the changes in the state matrix, A 0 , which reflects the new 
pole locations (as defined by the optimization), an eigenvalue assignment technique is 
employed. Specifically, the eigenvalue assignment technique discussed in Ref [101 is 

used, which is a sequential algorithm well-suited for partial assignment of eigenvalues 
in large-order systems. 
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Assume that the optimizer requires a subset of system poles to be reassigned to 
\ d . Now, consider a system (A 0 , B), where B is an n x r random matrix representing 
an arbitrary input influence matrix. The matrix B can have as many columns one 
would like, however, a reasonable choice would be max(m , r ). The change in the state 
matrix that would reassign the poles to A d is given by the gain matrix F such that the 
eigenvalues of A 0 -BF arc assigned to A d . Since the matrix B has typically more than 
one column, the gain matrix F is not unique, i.e., there is freedom beyond eigenvalue 
assignment. This freedom is used to minimize the norm of the gam matrix, and 
hence, the effective changes in the state matrix. The reader is referred to Ref. [10] for 
a detailed description of the eigenvalue assignment algorithm. Once the gain matrix 
F is computed, the new system realization is given by (A 0 - BF , B, C, D), i.e., 

G(z k ) = C(z k I n -A 0 + BF)- 1 B + D (26) 

In the next section, an optimization problem is posed for this approach. 

3.1.1 Absolute Error Cost Function 

Using a measure of the absolute error as the cost function, the optimization problem 
for system identification enhancement is defined as follows: 
minimize J : 

J = ||G(zfc) - G(z k )\\ F ( 27 ) 

over 

A* 

subject to 

|A.| < 1 

Here, A., represents a subset of the poles of the identified system G(z k ) which 
are used as design parameters, and | | F denotes the Frobenius norm. The constraint 
on the modulus of A s guarantees the stability of the identified system, and can be 
omitted if stability is not of concern. The computational procedure is as follows: 

1. Compute new pole locations, A, using optimizer. 
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2. Determine the new state matrix, with eigenvalues at A s , 
eigenvalue assignment technique described in Ref. [10]. 


using the sequential 


3. Update the identified transfer function matrix G(z k ) from Eq. (26), and compute 
the cost function J from Eq. (27). 


e number of poles that can be used as design parameters in the optimization is 
arbitrary. One can use all the poles in the system, or just a few, for example, the real 
poles of the system. If one starts the optimization with a system realization from the 
least-squares solution of Eq. (13), then it is very likely that the eomplex poles of the 
.dent, lied system, representing resonant peaks in the frequency response plots, match 
he experimental results well, and hence need not be manipulated any further. In 
such a case, real poles of the system and unstable poles, real or complex, are the best 
candidates for design optimization. However, one could conceivably use the modulus 

o complex poles, which determine the damping associated with each mode as 
design parameters as well. 

One of the problems with nonlinear programming is the tendency of the solution to 
converge to a local minimum. The problem becomes more aggravated as the number 
of design parameters increases. One way to deal with this problem is to restart the 
optimization from another set of design points in the neighborhood of the last optimal 
esign. Another way of avoiding this problem is to introduce an additional constraint 
requiring that the cost function be less than a desired value, i.e., 


- d (28) 

This constraint would move many of the local minima to the infeasible region, thereby 
avoiding them. An additional technique for dealing with local minima, which is more 

Tw rl the ° Ptimi2ation probiem that is be “S solved, is to restart the optimization 
of Eq. (27) with a new fictitious input influence matrix. Let B, represent the first 

input influence matrix, then choose the second influence matrix B 2 randomly from 
the null space of B\, i.e., 

= ( 29 ) 
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where N T - B\ = 0, and Ry is a random matrix of appropriate dimensions. The transfer 
B\ 

function of the system is then computed as 

G(z k ) = C(z k I n - Ai + B 2 F 2 )~ l B + D (30) 

with 

A\ — Ao — B\F\ 

Now, the optimization of Eq. (27) is restarted with B 2 as the input influence matrix. 
This process can be continued until either an acceptable cost function is achieved 
or the computational burden of optimization outweighs any reductions in the cost 
function. 

The cost function in Eq. (27), which is the Frobcnius norm of the error in the 
transfer functions (experimental and identified), is dominated by the peaks (reso- 
nants) of the transfer functions. Hence, optimization with Eq. (27) works well m 
reducing the errors at or around those peaks, or wherever the transfer function mag- 
nitude is significant, but it may not do much in reducing the errors elsewhere, e.g., 
zeros. In fact, the errors around the valleys might become worst. An approach to 
deal with this problem is presented next. 

3.1.2 Relative Error Cost Function 

A more equitable trade between the errors for peaks and valleys can be obtained by 
considering a complementary optimization problem, wherein a norm of the relative 
error is optimized instead of the absolute error given in Eq. (27). The optimization 
problem is posed as 
minimize J 2 : 

J 2 = \\[G(z k )-G(z k )]./G(z k )\\ F (31) 

over 


subject to 


J < Jd 
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I'M < 1 

Here, J is the Frobenius norm of the absolute error given by Eq. (27) and ” /” 
denotes element by element division. The upper bound limit J d can be any desired 
value, but might be set at the optimal value of J, using a previously computed value 
from Eq. (27), or might be set at the value of J from an initial least-squares solution. 
Similarly, the initial values of the design parameters may be set to previously obtained 
optimal values from the optimization problem in Eq. (27), or may be set to the values 
obtained from a least-squares solution. 

The second nonlinear programming approach for least-squares solution improve- 
ment is described next. 

3.2 Parameter optimization: direct approach 

A subset of system poles, as well as some additional coefficients to adjust the columns 
(rows) of the input (output) influence matrix, is used as optimization parameters. 
The optimizer adjusts these parameters to minimize the difference between the ex- 
perimental transfer function and the identified transfer function over frequency range 
of interest. In this approach, the identified state matrix A is first transformed to 
a real Schur canonical form or modal form. The system poles, which reside on the 
diagonal elements or 2 x 2 block diagonal partitions of the state matrix, are directly 
changed via optimization. The additional design parameters include coefficients for 
a set of basis vectors, which affect the input (output) influence matrix. To minimize 
the number of optimization parameters, the smaller of the input or output influence 
matrices is chosen for parameterization. Assume that the input influence matrix is 
the smaller of the two, i.e., there are more outputs than inputs, and let (A, B 0 , C, D) 

represent an initial realization for the identified system. Then, parameterize the input 
influence matrix B as 

H = Ho S + Nb 0 R ^22) 

where the matrix N Bo is an n x q matrix whose column span the null space of the 
matrix H 0 , that is 

Ab 0 H o = 0 
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and S and R are r x r and q x r coefficient matrices, respectively. An optimiza- 
tion problem, using an absolute error cost function, is discussed next using these 
parameters. 

3.2.1 Absolute Error Cost Function: Direct Approach 

The optimization problem is similar to the previous formulation derived for the eigen- 
value assignment technique, and is summarized as 
minimize J$: 

J 3 = \\G{z k ) - G(z k )\\ F ( 33 ) 


over 

blkdiag(A), S, R 


subject to 

|A(blkdiag(A)| < 1 

Here, G{z k ) represents the system realization given by 


G(z k ) = C{z k I n -A)- 1 B + D 


(34) 


The term blkdiag( ) denotes a subset of the 2x2 block diagonal of ( ). The constraint 
on the modulus of A(blkdiag(A)) guarantees the stability of the identified system, and 
can be omitted if stability is not of concern. The number of poles (eigenvalues), that 
are used as design parameters, is specified the designer. One can use all the poles in 
the system, or just a few, for example, the real poles of the system. In real Schur 
form or modal form, the real poles (eigenvalues) are on the diagonal, while the pairs 
of complex conjugate poles (eigenvalues) reside in the 2 x 2 block diagonals, in the 


following form 


real(Xi) imag(Xi)' 
—imag(Xi) real(Xi) . 


(35) 


where the real( ) and imag{ ) denote the real and imaginary parts, respectively. Care 
must be taken during parameterization to ensure those elements of the 2 x 2 block 
diagonals are parameterized properly so that complex conjugate pair remains intact 


throughout the optimization process. 
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In the next section, a complementary optimization problem is described using a 
relative error cost function. 

3.2.2 Relative Error Cost Function: Direct Approach 

A complementary optimization problem, wherein a norm of the relative error is opti- 
mized, may be used in order to provide a more equitable trade between errors in the 

peak and valley regions of the transfer function matrix. This optimization is defined 
as 

Ja = \\[G(z k ) - G{z k )]./G{z k )\\ F (36) 

over 

blkdiag(A), S, R 

subject to 

J3 < Ji d 

|A(blkdiag(A)| < 1 

Here, J 3 is the Frobemus norm of the absolute error given by Eq. (33), J 3rf can be any 
esired value, but it is often set to a previously computed value for J 3 in Eq. (33), 
or the value of J from in initial least-squares solution. Similarly, the initial values 
of the design parameters may be set to the previously computed optimal values for 

the optimization problem of Eq. (33), or may be set to the values obtained from a 
least-squares solution. 

The third nonlinear programming approach for least-squares solution improvement 
is described next. 


3.3 Parameter optimization: Least-Squares 

This method starts with selecting a subset of system poles as optimization parameters 
to minimize the error between the experimental and the identified transfer functions 
over frequency range of interest. The optimizer reassigns the system poles, which 
reside on the diagonal elements or 2x2 block diagonal partitions of the state matrix. 
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At each optimization step, corrections are made to the B : C , and D matrices through 
two least-squaxes solutions. 

similar to the previous two nonlinear programming approaches, two optimization 
problems are presented in the next two sections, one using an absolute error cost 
function and the other using a relative error cost function. 


3.3.1 Absolute Error Cost Function: Least-Squares 

let (A, Bo, Co, D 0 ) represent an initial realization for the identified system, and pa- 
rameterize the input and output influence matrices as follows, 

D = B 0 Sb + N Bo Rd = [B 0 N„ | ( ) = BQg (37) 

C=S c Co + RcNc a = f,Sc r c)[n°} sQcV <38) 

where the columns of N Bo represent a set of basis vectors in the null space of B 0 , the 
rows of N Co represent a set of basis vectors in the null space of C 0 , and Q B , defined 
in terms of S B and R B , and Q c , defined in terms of S c and R c , arc the appropriate 
coefficient matrices. These coefficients are determined, at each optimization step, by 
solving least-squares-based corrections of the absolute error norm. The optimization 
problem is given as 
minimize J§: 

J 5 = \\G(z k ) - G(z k )\\ F (39) 


over 

blkdiag(A) 


subject to 


|A(blkdiag(A))| < 1 


The complex matrix G(z k ) represents a system realization given by 


G{z k )=C{z k I n -A)- l B + D 


(40) 


The constraint on the modulus of A(blkdiag(A)) guarantees the stability of the iden- 
tified system, and can be omitted if stability is not of concern. At each optimization 
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step, as a new state matrix A is defined, corrections are performed to the B and D 

matrices via a least-squares solution, followed by corrections to the C and D matrices. 
These solutions are defined next. 

First^let G(z k ) = G{z k ) - D, and repartition the n d x (m x r) transfer function 
matrix, G(z k ), columnwise, such that each column of the repartitioned (n d x m) x r 
matrix, G^, is associated with an input. 

Define G{z k ) = C 0 (z k I n - A)~ l B , repartition G{z k ) similar to G col to obtain 
and define the absolute error function as 


e Gcol GcoiQb 

Now, solve for Q B , in a least-squares sense, to obtain 

Qb = ( G col G ml ) ~ 1 G col G col 
Once, Qb is computed, then D is computed as 


( 41 ) 


( 42 ) 


D = n ( G(z k ) — Co(z k I n — A)~ 1 B j 
where //( ) denotes the mean over frequency points. 

To compute Q c , first define G(z k ) = G(z k ) - D, and repartition the n d x (m x r) 
transfer functionmatrix, G(z k ), rowwise, such that each row of the repartitioned m x 
(n d x r) matrix, G row , is associated with an output. Define G(z k ) = C(z k I n - A )' 1 B, 
repartition G(^ k ) similar to G roU) to obtain G row , and define the error function as 

e — Grow QcG row 

Now, solve for Q c , in a least-squares sense, to obtain 

= Grow G row (G row G row ) 1 
Once, Q c is computed, then D is recomputed as 

D = fi (G( 2fe ) - C(z k I n - A)~ l B) 

The number of poles (eigenvalues) that can be used as design parameters is somewhat 

arbitrary. One may use all the poles in the system, or just a few, for example, the 
real poles of the system. 


( 44 ) 

( 45 ) 
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3.3.2 Relative Error Cost Function: Least-Squares 

Similar to the previous cases, the cost function is written in terms of a norm of the 
relative error to provide a more equitable trade between errors in the peak and valley 
regions of the transfer function matrix. This optimization problem is defined as: 
minimize J$: 

j 6 = ||[g( 2 o-gM-/g(zOIIf < 47 ) 

over 

blkdiag(;4) 

subject to 

J<> < -h d 

|A(blkdiag(A)| < 1 

Here, J 5 is the Frobenius norm of the absolute error given by Eq. (39), J 5d can be any 
desired value, but it is often set to a previously computed value for J 5 in Eq. (39), 
or the value of J from in initial least-squares solution. Similarly, the initial values 
of the design parameters may be set to the previously computed optimal values for 
the optimization problem of Eq. (39), or may be set to t he values obtained from a 
least-squares solution. The transfer function matrix G(z k ) is the same as that used 
for the pervious optimization problem, and is provided by Eq. (40), and the input 
and output influence matrices B and C are expanded following Eqs. (37)-(38). 

Now, there are two possible approaches for performing the least-squares-based 
corrections to the B , C, and D matrices. The first approach is to use the proce- 
dure outlined earlier for computing the coefficient matrices Q B (Eq. (42)) and Q c 
(Eq. (45)), exactly. That procedure is based on a least-squares correction of the abso- 
lute error norm. The second approach is to compute these coefficient matrices from a 
least-squares correction of the relative error norm. However, the least-squares prob- 
lems, whose solutions yield the coefficient vectors Q B and Qc, are somewhat different 
because they use a relative error norm instead of an absolute error norm. Without 
going into details of this alternative, the equations arc fairly similar to those obtained 
for the absolute error case, with three exceptions. First, instead of two least-squares 
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problems, one for input influence matrix correction and one for output influence ma- 
trix correction, there would be r least-squares problems for the input influence matrix 
correction, one for each column of the B matrix, and m least-squares problems for 
the output influence matrix correction, one for each row of the C matrix. The second 
exception is that one column (row) of the matrix G col (G row ), at a time, is used for 
each least-squares problem, with all its elements replaced by ones. The last exception 
is that the columns (rows) of G^Gr^) are normalized by the elements of G ml (G row ). 

4 Applications 

This section describes the application of the proposed techniques to the system iden- 
tification for the PARTI wmd-tunnel model [12], a laboratory test structure at NASA 
Langley. The model is a five-foot long, high aspect ratio wing designed to flutter at 
low speeds to simplify aerodynamic analyses and wind-tunnel testing. The fully as- 
sembled semi-span model is shown in Fig. 1. The model has a total of 72 actuators 
bonded to both sides of the plate. Each actuator contains two stacks of two 0.01 
inch piezoelectric patches. The 72 actuators are hardwired to actuate in 15 different 
groups. The 15 groupings were chosen such that each group primarily affects one of 
the first three natural modes. Each group can be considered as one input, because all 
the actuators in the group use the same signal. The piezoelectric patches were only 
used for actuation; ten strain gages and four accelerometers were used as sensors. As 
a result, there are a total of 15 inputs and 14 outputs. 

4.1 Weighted Least-Squares Solution 

In the first application, the transfer function from input No. 1 to all outputs is 
considered for identification. With signals from 14 sensor outputs ( m = 14), input 
No. 1 (r = 1), and 10th order polynomials ( p = 10) used in the matrix-fraction 
expansion (see Eqs. (l)-(3)), a weighted least-squares solution was first obtained 
from Eq. (22), using an exponential weighting function, given as 

w k \ k = 0 , 1 ,... J 
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Here, k = 0 refers to the zero frequency component of the FRF often known as 
the direct current (DC) term in electrical engineering, and w was chosen at 0.98. 
By adjusting the value of w one may emphasize the low frequencies or the higher 
frequencies. Values of w less than 1 would emphasize the lower end of the frequency 
spectrum. Here, w was set to 0.98, to emphasize the range of frequency from 0 to 
25 Hz. The weighted least-squares solution resulted in an identified model of order 
140, which included 4 unstable poles. However, since the actual testbed is stable, it 
is desired to obtain a stable identified model. Truncating the unstable states yielded 
a 136-order state-space realization of the system. Magnitude and phase FRF plots 
for outputs No. 7 and 10 are shown in Figs. 2 and 3, respectively. 

Comparison of the plots indicates an excellent agreement between the experimen- 
tal FRF (solid line) and the identified FRF (dashed line), particularly around the 
peaks of the FRF or where the gain values are significant. However, discrepancies 
can be observed around some of the zeros as well as where the gain values are small. 
This should be expected as the least-squares problem is dominated by the peaks and 
large gain values. Further inspection of these plots also indicates that the agreement 
between the experimental and identified results is better in the 0-25 Hz range. The 
Frobenius norm of the error between the experimental and identified transfer func- 
tions was computed at 90.128, the majority of which is due to the differences between 
two FRFs at DC frequency. In fact, since the DC gain values are quite large, par- 
ticularly in some output channels, they tend to dominate the rest of the FRF in a 
least-squares solution. Keep in mind that the DC gain values may not be accurate 
due to the use of accelerometers and their insensitivity at very low frequencies, drift 
problems that hampers accurate measurements, and lack of sufficient data. There- 
fore, in this case, it is reasonable to de-emphasize the DC values by assigning a zero 
value to the corresponding weighting function, such that the DC weight is set to 
zero. The Bode plots, using the modified weighting function , are shown in Figs. 2 and 
3 as dashed-dotted lines. These figures indicate moderate improvements in various 
frequency ranges. The Frobenius norm of the error between the experimental and 
identified transfer functions was computed at 90.134, a very minor change from the 
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previous results, indicating that little in terms of minimizing the overall error is lost. 
Comparing, the norm of the error for all frequency points except DC, shows that the 
error went down from 0.241 to 0.223, which quantifies the better match by using the 
modified weighting function. 

In order to show the effectiveness of the modified exponential weighting, a polyno- 
mial with p = 3 is used in the matrix-fraction description. The weighted least-squares 
solution resulted in a stable identified model of order 40. Figures 4 and 5 illustrate 
the stable least-squares solutions for the nominal and modified exponential weight- 
ings for outputs No. 7 and 10. These figures clearly demonstrate the advantage of 
modified exponential weighting for this problem. In fact, the Frobenius norm of the 
error between the experimental and identified transfer functions dropped from 12.035 
to 0.3021, a significant improvement. 

4.2 Further Enhancements: Nonlinear Programming 

To demonstrate the potential of the nonlinear programming approaches to further 
enhance the least-squares solution, the direct optimization approach is applied to an 
identified model for the PARTI testbed, obtained from a least-squares solution. Three 
optimization problems were considered and carried out in this application. The first 
optimization approach was the one posed in Eq. (33), minimizing the Frobenius norm 
of the absolute error cost function for all frequency points except the DC, using as 
design variables the real poles of the identified plant as well as the coefficients of the 
basis vectors, 5 and R, from Eq. (32). The initial least-squares solution provided a 
stable model of order 40, with 6 real poles. This model was transformed into modal 
form to make it amenable to optimization with the direct approach. There were 26 
design variables used in the optimization, the first 6 being the values of the real poles 
of the system, the 7th representing the coefficient S (see Eq. (32)), and the remaining 
19 representing the elements of the coefficient vector R (see Eq. (32)), corresponding 
to 19 basis vectors from the null space of the matrix B 0 . It should be noted that one 
could have included the complex poles of the system as design variables too. However, 
it was concluded that this may not be necessary since the dominant poles (resonances) 
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of the system were captured by the least-squares solution, as observed from Figs. 4 
and 5. The optimization included 7 constraints, the first six to assure the stability 
of the system as the poles were reassigned, and the last constraint on the value of 
the error norm to avoid convergence to undesirable local minima. The optimization 
was performed using the Automated Design Synthesis (ADS) software [13]. The 
interior penalty function method of ADS was used to solve the nonlinear programming 
problems. In this method, the constrained optimization problem is transformed into 
an unconstrained problem through creation of a objective function which is the sum of 
the original objective function and an imposed penalty function (which is a function 
of the constraints) [14]. The initial solution used in the optimization was that of 
stable least-squares solution with modified exponential weighting. The optimization 
reduced the norm of the absolute error from the initial value of 0.302 to 0.285. The 
optimization was repeated once with a new input influence matrix bases. These bases 
included the initial influence vector, along with a new set of 19 random basis vectors 
from the null space of the matrix [Bq Nb 0 } ■ The second optimization run reduced 
the norm of the error from 0.289 to 0.249, which is fairly close to 0.223, the error 
norm corresponding to the stable least-squares 'solution with p = 10. This process 
could have been continued further, however, it was deemed that the computational 
burden would outweigh any additional reductions in the norm of the error. The Bode 
plots, comparing the transfer function matrices for the experimental, nominal, and 
optimal results are provided in Figs. 6 and 7, for output No. 7 and 10, where it 
is observed that the identified model, via optimization, performs well, almost at a 
comparable level to that for p = 10. This is very encouraging, considering that the 
optimization-based model is a 40th order model compared to 136th order model for 

p - 10. 

Figures 6 and 7 also indicate that there is room for improving the test and analysis 
agreement around the transmission zeros and low gain regions. To accomplish this, a 
second optimization problem was performed. This case used the formulation described 
in Eq. (36), in which the Frobenius norm of the relative error is minimized. The initial 
design used the optimal solution obtained from the first optimization problem. Similar 
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to the first optimization, 26 design variables were used in the optimization, the first 
6 being the location of the real poles of the system; the 7th represents the coefficient 
S (see Eq. (32)); and the remaining 19 representing the elements of the coefficient 
vector R (see Eq. (32)), corresponding to 19 basis vectors from the null space of the 
matrix B 0 . In addition, 7 constraints were used, the first six to guarantee system 
stability, and the last constraint to ensure that the Probenius norm of the absolute 
error remained equal to or below 0.249, the level obtained by the first optimization. 
The optimization reduced the norm of the relative error from the initial value of 
181.11 to 152.48. The optimization was repeated once more with a set of new input 
influence matrix bases. The second optimization run reduced the norm of the error 
from 152.48 to 118.34, and the optimization process was terminated. The Bode plots 
for the experimental, nominal, and optimal results, for outputs No. 7 and 10, are 
presented in Figs. 8 and 9. These figures illustrates a considerable improvement in 
the optimal model in matching the area around the transmission zeros and low gain 
regions of the transfer functions. These results demonstrate the advantages of further 
enhancements to the least-squares solution. 

Another optimization problem considered was the problem posed in Eq. (39) where 
the Frobenius norm of the absolute error is minimized while adjusting the eigenvalues 
of the state matrix, subject to stability constraints. Moreover, the optimization in- 
cluded corrections to the B and D matrices, followed by corrections to the C and D 
matrices, at each functional evaluation (see Eq. (42)-(43) and (45)-(46)). The 6 de- 
sign variables used in the optimization were the values of the real poles of the system. 
The optimization included 7 constraints, the first six to guarantee the stability of the 
systems as the poles were reassigned, and a constraint on the value of the error norm 
to avoid undesirable local minima. The initial design used in the optimization was 
was taken from a stable least-squares solution with modified exponential weighting, 
i.e., zero DC weighting. The optimization reduced the norm of the absolute error 
from the initial value of 0.250 to 0.197, which is over 20% reduction. Plots, compar- 
ing the transfer function matrices for the experimental, nominal, and optimal results 
are provided in Figs. 10 and 11, for output No. 7 and 11. The identified model (via 
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optimization) agrees very well with the experimental data. Also note that output 1 1 
is now shown as opposed to 10 in the previous cases. This output sensor shows modes 
around 28 and 38 Hz which are also being identified. 

All the identification results obtained so fax were based on the 14 transfer functions 
from the first input to all 14 outputs. Now, let us consider the transfer functions from 
all 15 inputs to all 14 outputs for identification. With the signals from all 14 sensor 
outputs (tyi = 14) and all 15 inputs (r = 15), and 3rd order polynomials ( p 
3) used in the matrix-fraction expansion(see Eqs. (l)-(3)), a weighted least square 
solution was first obtained from Eq. (22). Similar to the previous cases, an exponential 
weighting function was used, with parameter w chosen at 0.98 to emphasize the range 
of low frequencies. In addition, the DC weight was set to zero. The Probenius norm 
of the error between the experimental and identified transfer functions was computed 
at 246.855, the majority of which is due to a discrepancy between two FRFs at the 
DC frequency, i.e., zero frequency. The Frobenius norm of the absolute error, for 
all frequency points except the DC, was computed to be 1.721. For the purpose of 
illustration, plots for the experimental and realized transfer functions are depicted in 
Figs. 12 and 13, for outputs No. 7 and 11 with input No. 1, and in Figs. 14 and 
15, for the same outputs with input No. 8. The experimental transfer functions arc 
shown as solid lines and the transfer functions, obtained via direct least square, as 
dashed lines. These figures indicate moderate to good agreement between the transfer 
functions in low frequencies ranges, particularly, around the peaks or high gain areas 
of the transfer functions. 

Now consider the least-squares optimization approach presented in Eq. (39) for 
the 15 inputs and 14 output case. The initial design used in the optimization was the 
stable least-squares solution with modified exponential weighting. This realization 
had 14 real poles, whose locations were used as design variables in the optimization, 
i.e., there were 14 design variables. The optimization included 15 constraints, the first 
14 to guarantee the stability of the systems as the poles were reassigned, and the last 
constraint on the value of the error norm to avoid convergence to undesirable local 
minima. The optimization reduced the norm of the error from the initial value of 
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1.165 to 1.090, about a 6.5% reduction. Plots, comparing for the experimental (solid 
line), nominal (dashed line), and optimal transfer function (dashed-dotted line) are 
provided in Figs. 12 through 15. It is observed that the identified model (obtained 
via optimization) performs well, although only 3rd order polynomials was used in the 
matrix fraction description to match a 15 input by 14 output transfer function. Com- 
parison of Figs. 12 and 13 with Figs. 10 and 11, which correspond to the same input 
and output channels, and were obtained from a single input identification problem, 

reconfirms the good level of correlation obtained following the optimization-based 
approach. 

5 Concluding Remarks 

Several techniques have been presented to identify an experimental system model 
directly from frequency response data. The techniques used a matrix-fraction de- 
scription (MFD) to describe the identified system. The MFD coefficients were ob- 
tained from the solution of a weighted least-squares problem. Frequency weighting 
concepts were introduced in order to emphasize a frequency range of interest. Three 
optimization-based methods were introduced to fine-tune the experimentally realized 
models. The first method uses an eigenvalue assignment technique to reassign a sub- 
set of the system poles to enhance the approximation of the identified model. The 
second method considers the model in the real Schur or modal form, and uses a subset 
of system poles as well as additional parameters to adjust the columns (rows) of the 
input (output) influence matrix in order to improve the model fit. The third method 
optimizes a subset of the system poles, while the input and output influence matri- 
ces are computed at every optimization step through least-squares procedure. The 
methods were applied to data from PARTI wind tunnel model, a laboratory testbed 
at NASA Langley Research Center. The benefits of the optimization-based refine- 
ment techniques as well as frequency weighting techniques were demonstrated. It was 
shown that with optimal fine-tuning and proper choice of frequency weighting a 40th 
order system realization could provide almost the same level of model fit as a full-order 
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136th order model. Specifically, the third optimization method provided the largest 
improvement in the fit of the experimental data, with a 20-perccnt smaller absolute 
error norm than that obtained from the second optimization method. However, the 
computational cost of performing the third optimization formulation was considerably 
higher than the second optimization problem, primarily due to the large least-squares 
solutions that were solved at every optimization step to perform corrections to the 
B, C , and D matrices. The optimizations used finite difference techniques to provide 
gradient information on the objective function and constraints. The numerical com- 
putation of the gradients may require a large number of functional evaluation, which 
would be costly in a computational sense. Alternatively, one may attempt to obtain 
analytical expressions for the gradients, and perhaps second-order partial derivatives, 
to improve computational efficiency and accuracy. 
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Figure 1: PARTI Model 
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Output No. 7 




Figure 2: Comparison of FRFs for Output No. 7 with exponential weighting and 136- 

FRFwRh C n P r rimen u tal FI } F (S ° lid llne)j identified FRF (dashed line), identified 
FRF with zero DC weighting (dashed- dotted line). 
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Output No. 10 



Frequency (Hz) 



Frequency (Hz) 

Figure 3: FRF for Output No. 10 with exponential weighting and 136-order system: 
experimental FRF (solid line), identified FRF (dashed line), identified FRF with zero 
DC weighting (dashed-dotted line) 
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Output No. 7 




Figure 4: Comparison of FRFs for Output No. 7 with exponential weighting and 
40-order of system: experimental FRF (solid line), identified FRF (dashed line), 
identified FRF with zero DC weighting (dashed-dotted line) 
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Output No. 10 




Figure 5: Comparison of FRFs for Output No. 10 with exponential weighting and 
40-order of system: experimental FRF (solid line), identified FRF (dashed line), 
identified FRF with zero DC weighting (dashed-dotted line) 
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Fi ^ 6: Comparison of FRFs for Output No. 7 with direct optimization approach 
and 40-order system: experimental FRF (solid line), identified FRF (dashed line) with 
zero DC weighting, enhanced FRF with absolute-error optimization (dashed-dotted 
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Phase (deg) Amplitude 



Figure 7: Comparison of FRFs for Output No. 10 with direct optimization approach 
and 40-order system: experimental FRF (solid line), identified FRF (dashed line) with 
zero DC weighting, enhanced FRF with absolute-error optimization (dashed-dotted 

line) 
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and 40-order system: experimental FRF (solid line), identified FRF Mash ’d lj„ e , 
with zero DC weighting; enhanrpr! prp 1 v Q iinej 

dotted line) enhanced FRF with relative-error optimization (dashed- 
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with zero DC weighting, enhanced FRF with relative-error optimization (dashed- 
dotted line) 
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Figure 10: Comparison of FRFs for Output No. 7 with least-squares optimiza- 
tion approach and 40-order system: experimental FRF (solid line), identified FRF 
(dashed line) with zero DC weighting, enhanced FRF with absolute-error optimiza- 
tion (dashed-dotted line) 
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Figure 11: Comparison of FRFs for Output No. 11 with least-squares optimiza- 
tion approach and 40-order system: experimental FRF (solid line), identified FRF 
(dashed line) with zero DC weighting, enhanced FRF with absolute-error optimiza- 
tion (dashed-dotted line) 
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Figure 12: Comparison of FRFs for Output No. 7 and Input No. 1 with least-square 

& u P JT Ch 3X1(1 42 ' order s y stem: experimental FRF (solid line) identi 
fied FRF (dashed line) with zero DC weighting, enhanced FRF with absolute-erro 
optimization (dashed-dotted line) 


40 





Phase (deg) Amplitude 




Figure 13: Comparison of FRFs for Output No. 11 and Input No. 1 with least- 
squares optimization approach and 42-order system: experimental FRF (solid line), 
identified FRF (dashed line) with zero DC weighting, enhanced FRF with absolute- 
error optimization (dashed-dotted line) 
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Figure 14: Comparison of FRFs for Output No. 7 and Input No. 8 with least-squares 
fied FRF and 42 ‘° rder s y stem: experimental FRF (solid line) klenti- 
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Frequency (Hz) 


Figure 15: Comparison of FRFs for Output No. 11 and Input No. 8 with least- 
squares optimization approach and 42-order system: experimental FRF (solid line), 
identified FRF (dashed line) with zero DC weighting, enhanced FRF with absolute- 
error optimization (dashed-dotted line) 


43 




REPORT DOCUMENTATION PAGE 


Form Approved 
OMBNo. 0704-0188 


p w — I o'/vfo mo. U/U4—U188 

^thwinj and malmnnj the data naaded. ‘"‘S’ ’’’* far nKm, ' n 9 inductions, aoarchmg exist, nfl data oourcea 

r^'Tn MUamB SU0 « ,Mi " Ons *>' W. burden, 'to Washington H«dq«™ V"’ ^ W "" ° ,her “<*< <* »"« 

P T ouency USE mu v T^ZZ?**- - - ,he | ^'“rrmiMim 1 — Bu a9 "' "T ,2 ’ * Da ™ 


2. REPORT DATE 

April 1999 

| 4. TITLE AND SUBTITLE 

Optimal Frequency-Domain System Realization with Weighting 


1 3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 

5. FUNDING NUMBERS 

WU 522-32-41-03 


| 6 AUTHOR(S) 

Jer-Nan Juang and Peiman G. Maghami 


I 7 PERFORMING ORGANIZATION NAME(S) AND ADDRESSEES) 

NASA Langley Research Center 
Hampton, VA 23681-2199 


* 3tJ(JN5 OR'NG/MONITORINO AG ENCY NAME(S) AND ADDRESSES) 

I National Aeronautics and Space Administration 
Washington, DC 20546-0001 


1 11. SUPPLEMENTARY NOTES 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

L-17845 


m SPONSOHING/MONfTORING 
AGENCY REPORT NUMBER 

NASA/TM-1 999-2091 35 


12b. DISTRIBUTION CODE 


M2a. DISTRIBUTION/AVAI LABILITY STATEMENT — 

Unclassified-Unlimited 

Subject Category 39 Distribution: Standard 

Availability: NASA CASI (301) 621-0390 


1 13. ABSTRACT (Umximum 200 word*) 

^ *“ *•*■"* response 0s, a. 

weighting is introduced to solve a weighted least-squares DroW^T^t. Fr ^" C !7 e ' 9h “ n8 such as “‘Ponential 
matrix-fraction description. A multi-variable t °° bta,n the coefficient matrices for the 

the matrix-fraction description Three different aonmarhT* & ? formed usln 9 the coefficient matrices of 

programming methods to'SLSS.' "* ^ usinfl 
technique to reassign a subset of system dX mTmnrnl ?k Jh " ,eth ° d , uses an eigenvalue assignment 

model in the real Schur or modal form reassions a S uh«Jr If ld ® ntlf,ed model - The second method deals with the 
input (output) influence matrix using a nonlinearoptfmfzer The^frd meh^T columns < rows > of the 

the input and output influence mathces are rst'n^ar^^^I^mlMfion's^p tmo^h te'^h^uamspro^clures 111 


14. subject tpquc 

System Identification, System Realization, Modal Parameter Identification 
Frequency Response Function 

Ts. NUMBER OF PAGES 

48 

17. SECURmr CLASSIRCATION 
OF REPORT 

Unclassified 
NSN 7540-01-280-5500 

18. SECURITY CLASSIRCATION 
OF THIS PAGE 

Unclassified 

19. SECURITY CLASSIRCATION 
OF ABSTRACT 

Unclassified 

16. PRICE CODE 

A04 

20. LIMITATION OF ABSTRACT 


Prescribed by ANSI Sid. Z39-18 
298-102 


